From d199e0010a66d63ed407a1f194f73dad078840b3 Mon Sep 17 00:00:00 2001 From: Ell Date: Sun, 10 Sep 2017 11:30:56 -0400 Subject: [PATCH] babl: add babl-polynomial.[hc] Implements BablPolynomial, an opaque type representing a real valued polynomial of a real variable. Currently, exposes the following functions: * babl_polynomial_eval(): Evaluates a polynomial. * babl_polynomial_approximate_gamma(): Calculates a polynomial approximation of a gamma function. --- babl/Makefile.am | 2 + babl/babl-internal.h | 1 + babl/babl-polynomial.c | 529 +++++++++++++++++++++++++++++++++++++++++ babl/babl-polynomial.h | 79 ++++++ 4 files changed, 611 insertions(+) create mode 100644 babl/babl-polynomial.c create mode 100644 babl/babl-polynomial.h diff --git a/babl/Makefile.am b/babl/Makefile.am index 9fddfbf..bc13d95 100644 --- a/babl/Makefile.am +++ b/babl/Makefile.am @@ -29,6 +29,7 @@ c_sources = \ babl-model.c \ babl-mutex.c \ babl-palette.c \ + babl-polynomial.c \ babl-ref-pixels.c \ babl-sampling.c \ babl-sanity.c \ @@ -61,6 +62,7 @@ h_sources = \ babl-model.h \ babl-matrix.h \ babl-mutex.h \ + babl-polynomial.h \ babl-ref-pixels.h \ babl-sampling.h \ babl-space.h \ diff --git a/babl/babl-internal.h b/babl/babl-internal.h index fcdfe16..7bf0544 100644 --- a/babl/babl-internal.h +++ b/babl/babl-internal.h @@ -51,6 +51,7 @@ #include "babl-memory.h" #include "babl-mutex.h" #include "babl-cpuaccel.h" +#include "babl-polynomial.h" /* fallback to floor function when rint is not around */ #ifndef HAVE_RINT diff --git a/babl/babl-polynomial.c b/babl/babl-polynomial.c new file mode 100644 index 0000000..d51685d --- /dev/null +++ b/babl/babl-polynomial.c @@ -0,0 +1,529 @@ +/* babl - dynamically extendable universal pixel conversion library. + * Copyright (C) 2017, Øyvind Kolås and others. + * + * babl-polynomial.c + * Copyright (C) 2017 Ell + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 3 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General + * Public License along with this library; if not, see + * . + */ + +#ifdef BABL_POLYNOMIAL_DEGREE + +BABL_POLYNOMIAL_DEGREE ( 0, __) +BABL_POLYNOMIAL_DEGREE ( 1, 0) +BABL_POLYNOMIAL_DEGREE ( 2, 1) +BABL_POLYNOMIAL_DEGREE ( 3, 2) +BABL_POLYNOMIAL_DEGREE ( 4, 3) +BABL_POLYNOMIAL_DEGREE ( 5, 4) +BABL_POLYNOMIAL_DEGREE ( 6, 5) +BABL_POLYNOMIAL_DEGREE ( 7, 6) +BABL_POLYNOMIAL_DEGREE ( 8, 7) +BABL_POLYNOMIAL_DEGREE ( 9, 8) +BABL_POLYNOMIAL_DEGREE (10, 9) +BABL_POLYNOMIAL_DEGREE (11, 10) +BABL_POLYNOMIAL_DEGREE (12, 11) +BABL_POLYNOMIAL_DEGREE (13, 12) +BABL_POLYNOMIAL_DEGREE (14, 13) +BABL_POLYNOMIAL_DEGREE (15, 14) +BABL_POLYNOMIAL_DEGREE (16, 15) +BABL_POLYNOMIAL_DEGREE (17, 16) +BABL_POLYNOMIAL_DEGREE (18, 17) +BABL_POLYNOMIAL_DEGREE (19, 18) +BABL_POLYNOMIAL_DEGREE (20, 19) +BABL_POLYNOMIAL_DEGREE (21, 20) +BABL_POLYNOMIAL_DEGREE (22, 21) + +#undef BABL_POLYNOMIAL_DEGREE + +#else + +#include "config.h" +#include +#include +#include "babl-internal.h" + + +#define BABL_BIG_POLYNOMIAL_MAX_DEGREE (2 * BABL_POLYNOMIAL_MAX_DEGREE + BABL_POLYNOMIAL_MAX_SCALE) +#define EPSILON 1e-10 + + +typedef struct +{ + BablPolynomialEvalFunc eval; + int degree; + int scale; + double coeff[BABL_BIG_POLYNOMIAL_MAX_DEGREE + 1]; +} BablBigPolynomial; + + +#define BABL_POLYNOMIAL_EVAL___(poly, x) 0.0 +#define BABL_POLYNOMIAL_EVAL_0(poly, x) ( (poly)->coeff[0]) +#define BABL_POLYNOMIAL_EVAL_1(poly, x) ( (poly)->coeff[1]) +#define BABL_POLYNOMIAL_EVAL_2(poly, x) (BABL_POLYNOMIAL_EVAL_0 (poly, x) * (x) + (poly)->coeff[2]) +#define BABL_POLYNOMIAL_EVAL_3(poly, x) (BABL_POLYNOMIAL_EVAL_1 (poly, x) * (x) + (poly)->coeff[3]) +#define BABL_POLYNOMIAL_EVAL_4(poly, x) (BABL_POLYNOMIAL_EVAL_2 (poly, x) * (x) + (poly)->coeff[4]) +#define BABL_POLYNOMIAL_EVAL_5(poly, x) (BABL_POLYNOMIAL_EVAL_3 (poly, x) * (x) + (poly)->coeff[5]) +#define BABL_POLYNOMIAL_EVAL_6(poly, x) (BABL_POLYNOMIAL_EVAL_4 (poly, x) * (x) + (poly)->coeff[6]) +#define BABL_POLYNOMIAL_EVAL_7(poly, x) (BABL_POLYNOMIAL_EVAL_5 (poly, x) * (x) + (poly)->coeff[7]) +#define BABL_POLYNOMIAL_EVAL_8(poly, x) (BABL_POLYNOMIAL_EVAL_6 (poly, x) * (x) + (poly)->coeff[8]) +#define BABL_POLYNOMIAL_EVAL_9(poly, x) (BABL_POLYNOMIAL_EVAL_7 (poly, x) * (x) + (poly)->coeff[9]) +#define BABL_POLYNOMIAL_EVAL_10(poly, x) (BABL_POLYNOMIAL_EVAL_8 (poly, x) * (x) + (poly)->coeff[10]) +#define BABL_POLYNOMIAL_EVAL_11(poly, x) (BABL_POLYNOMIAL_EVAL_9 (poly, x) * (x) + (poly)->coeff[11]) +#define BABL_POLYNOMIAL_EVAL_12(poly, x) (BABL_POLYNOMIAL_EVAL_10 (poly, x) * (x) + (poly)->coeff[12]) +#define BABL_POLYNOMIAL_EVAL_13(poly, x) (BABL_POLYNOMIAL_EVAL_11 (poly, x) * (x) + (poly)->coeff[13]) +#define BABL_POLYNOMIAL_EVAL_14(poly, x) (BABL_POLYNOMIAL_EVAL_12 (poly, x) * (x) + (poly)->coeff[14]) +#define BABL_POLYNOMIAL_EVAL_15(poly, x) (BABL_POLYNOMIAL_EVAL_13 (poly, x) * (x) + (poly)->coeff[15]) +#define BABL_POLYNOMIAL_EVAL_16(poly, x) (BABL_POLYNOMIAL_EVAL_14 (poly, x) * (x) + (poly)->coeff[16]) +#define BABL_POLYNOMIAL_EVAL_17(poly, x) (BABL_POLYNOMIAL_EVAL_15 (poly, x) * (x) + (poly)->coeff[17]) +#define BABL_POLYNOMIAL_EVAL_18(poly, x) (BABL_POLYNOMIAL_EVAL_16 (poly, x) * (x) + (poly)->coeff[18]) +#define BABL_POLYNOMIAL_EVAL_19(poly, x) (BABL_POLYNOMIAL_EVAL_17 (poly, x) * (x) + (poly)->coeff[19]) +#define BABL_POLYNOMIAL_EVAL_20(poly, x) (BABL_POLYNOMIAL_EVAL_18 (poly, x) * (x) + (poly)->coeff[20]) +#define BABL_POLYNOMIAL_EVAL_21(poly, x) (BABL_POLYNOMIAL_EVAL_19 (poly, x) * (x) + (poly)->coeff[21]) +#define BABL_POLYNOMIAL_EVAL_22(poly, x) (BABL_POLYNOMIAL_EVAL_20 (poly, x) * (x) + (poly)->coeff[22]) + +#define BABL_POLYNOMIAL_DEGREE(i, i_1) \ + static double \ + babl_polynomial_eval_1_##i (const BablPolynomial *poly, \ + double x) \ + { \ + const double x2 = x * x; \ + (void) x2; \ + return BABL_POLYNOMIAL_EVAL_##i (poly, x2) + \ + BABL_POLYNOMIAL_EVAL_##i_1 (poly, x2) * x; \ + } +#include "babl-polynomial.c" + +#define BABL_POLYNOMIAL_DEGREE(i, i_1) \ + static double \ + babl_polynomial_eval_2_##i (const BablPolynomial *poly, \ + double x) \ + { \ + return BABL_POLYNOMIAL_EVAL_##i (poly, x) + \ + BABL_POLYNOMIAL_EVAL_##i_1 (poly, x) * sqrt (x); \ + } +#include "babl-polynomial.c" + +static const BablPolynomialEvalFunc babl_polynomial_eval_funcs[BABL_POLYNOMIAL_MAX_SCALE] + [BABL_BIG_POLYNOMIAL_MAX_DEGREE + 1] = +{ + { + #define BABL_POLYNOMIAL_DEGREE(i, i_1) babl_polynomial_eval_1_##i, + #include "babl-polynomial.c" + }, + { + #define BABL_POLYNOMIAL_DEGREE(i, i_1) babl_polynomial_eval_2_##i, + #include "babl-polynomial.c" + } +}; + + +static void +babl_polynomial_set_degree (BablPolynomial *poly, + int degree, + int scale) +{ + babl_assert (degree >= BABL_POLYNOMIAL_MIN_DEGREE && + degree <= BABL_BIG_POLYNOMIAL_MAX_DEGREE); + babl_assert (scale >= BABL_POLYNOMIAL_MIN_SCALE && + scale <= BABL_POLYNOMIAL_MAX_SCALE); + + poly->eval = babl_polynomial_eval_funcs[scale - 1][degree]; + poly->degree = degree; + poly->scale = scale; +} + +static double +babl_polynomial_get (const BablPolynomial *poly, + int i) + +{ + return poly->coeff[poly->degree - i]; +} + +static void +babl_polynomial_set (BablPolynomial *poly, + int i, + double c) + +{ + poly->coeff[poly->degree - i] = c; +} + +static void +babl_polynomial_copy (BablPolynomial *poly, + const BablPolynomial *rpoly) +{ + poly->eval = rpoly->eval; + poly->degree = rpoly->degree; + poly->scale = rpoly->scale; + memcpy (poly->coeff, rpoly->coeff, (rpoly->degree + 1) * sizeof (double)); +} + +static void +babl_polynomial_reset (BablPolynomial *poly, + int scale) +{ + babl_polynomial_set_degree (poly, 0, scale); + babl_polynomial_set (poly, 0, 0.0); +} + +static void +babl_polynomial_shrink (BablPolynomial *poly) +{ + int i; + + for (i = 0; i <= poly->degree; i++) + { + if (fabs (poly->coeff[i]) > EPSILON) + break; + } + + if (i > 0) + { + memmove (poly->coeff, &poly->coeff[i], + (poly->degree - i + 1) * sizeof (double)); + + babl_polynomial_set_degree (poly, poly->degree - i, poly->scale); + } +} + +static void +babl_polynomial_add (BablPolynomial *poly, + const BablPolynomial *rpoly) +{ + int i; + + babl_assert (poly->scale == rpoly->scale); + + if (poly->degree >= rpoly->degree) + { + for (i = 0; i <= rpoly->degree; i++) + { + babl_polynomial_set (poly, i, babl_polynomial_get (poly, i) + + babl_polynomial_get (rpoly, i)); + } + } + else + { + int orig_degree = poly->degree; + + babl_polynomial_set_degree (poly, rpoly->degree, poly->scale); + + for (i = 0; i <= orig_degree; i++) + { + babl_polynomial_set (poly, i, poly->coeff[orig_degree - i] + + babl_polynomial_get (rpoly, i)); + } + + for (; i <= rpoly->degree; i++) + babl_polynomial_set (poly, i, babl_polynomial_get (rpoly, i)); + } +} + +static void +babl_polynomial_sub (BablPolynomial *poly, + const BablPolynomial *rpoly) +{ + int i; + + babl_assert (poly->scale == rpoly->scale); + + if (poly->degree >= rpoly->degree) + { + for (i = 0; i <= rpoly->degree; i++) + { + babl_polynomial_set (poly, i, babl_polynomial_get (poly, i) - + babl_polynomial_get (rpoly, i)); + } + } + else + { + int orig_degree = poly->degree; + + babl_polynomial_set_degree (poly, rpoly->degree, poly->scale); + + for (i = 0; i <= orig_degree; i++) + { + babl_polynomial_set (poly, i, poly->coeff[orig_degree - i] - + babl_polynomial_get (rpoly, i)); + } + + for (; i <= rpoly->degree; i++) + babl_polynomial_set (poly, i, -babl_polynomial_get (rpoly, i)); + } +} + +static void +babl_polynomial_scalar_mul (BablPolynomial *poly, + double a) +{ + int i; + + for (i = 0; i <= poly->degree; i++) + poly->coeff[i] *= a; +} + +static void +babl_polynomial_scalar_div (BablPolynomial *poly, + double a) +{ + int i; + + for (i = 0; i <= poly->degree; i++) + poly->coeff[i] /= a; +} + +static void +babl_polynomial_mul_copy (BablPolynomial *poly, + const BablPolynomial *poly1, + const BablPolynomial *poly2) +{ + int i; + int j; + + babl_assert (poly1->scale == poly2->scale); + + babl_polynomial_set_degree (poly, poly1->degree + poly2->degree, poly1->scale); + + memset (poly->coeff, 0, (poly->degree + 1) * sizeof (double)); + + for (i = poly1->degree; i >= 0; i--) + { + for (j = poly2->degree; j >= 0; j--) + { + babl_polynomial_set (poly, i + j, babl_polynomial_get (poly, i + j) + + babl_polynomial_get (poly1, i) * + babl_polynomial_get (poly2, j)); + } + } +} + +static void +babl_polynomial_integrate (BablPolynomial *poly) +{ + int i; + + babl_polynomial_set_degree (poly, poly->degree + poly->scale, poly->scale); + + for (i = 0; i <= poly->degree - poly->scale; i++) + { + poly->coeff[i] *= poly->scale; + poly->coeff[i] /= poly->degree - i; + } + + for (; i <= poly->degree; i++) + poly->coeff[i] = 0.0; +} + +static void +babl_polynomial_gamma_integrate (BablPolynomial *poly, + double gamma) +{ + int i; + + babl_polynomial_set_degree (poly, poly->degree + poly->scale, poly->scale); + + gamma *= poly->scale; + + for (i = 0; i <= poly->degree - poly->scale; i++) + { + poly->coeff[i] *= poly->scale; + poly->coeff[i] /= poly->degree - i + gamma; + } + + for (; i <= poly->degree; i++) + poly->coeff[i] = 0.0; +} + +static double +babl_polynomial_inner_product (const BablPolynomial *poly1, + const BablPolynomial *poly2, + double x0, + double x1) +{ + BablBigPolynomial temp; + + babl_polynomial_mul_copy ((BablPolynomial *) &temp, poly1, poly2); + babl_polynomial_integrate ((BablPolynomial *) &temp); + + return babl_polynomial_eval ((BablPolynomial *) &temp, x1) - + babl_polynomial_eval ((BablPolynomial *) &temp, x0); +} + +static double +babl_polynomial_gamma_inner_product (const BablPolynomial *poly, + double gamma, + double x0, + double x1) +{ + BablBigPolynomial temp; + + babl_polynomial_copy ((BablPolynomial *) &temp, poly); + babl_polynomial_gamma_integrate ((BablPolynomial *) &temp, gamma); + + return babl_polynomial_eval ((BablPolynomial *) &temp, x1) * pow (x1, gamma) - + babl_polynomial_eval ((BablPolynomial *) &temp, x0) * pow (x0, gamma); +} + +static double +babl_polynomial_norm (const BablPolynomial *poly, + double x0, + double x1) +{ + return sqrt (babl_polynomial_inner_product (poly, poly, x0, x1)); +} + +static void +babl_polynomial_normalize (BablPolynomial *poly, + double x0, + double x1) +{ + double norm; + + norm = babl_polynomial_norm (poly, x0, x1); + + if (norm > EPSILON) + babl_polynomial_scalar_div (poly, norm); +} + +static void +babl_polynomial_project_copy (BablPolynomial *poly, + const BablPolynomial *rpoly, + const BablPolynomial *basis, + int basis_n, + double x0, + double x1) +{ + int i; + + babl_assert (basis_n > 0); + + babl_polynomial_reset (poly, basis[0].scale); + + for (i = 0; i < basis_n; i++) + { + BablPolynomial temp; + + babl_polynomial_copy (&temp, &basis[i]); + babl_polynomial_scalar_mul (&temp, + babl_polynomial_inner_product (&temp, rpoly, + x0, x1)); + babl_polynomial_add (poly, &temp); + } +} + +static void +babl_polynomial_gamma_project_copy (BablPolynomial *poly, + double gamma, + const BablPolynomial *basis, + int basis_n, + double x0, + double x1) +{ + int i; + + babl_assert (basis_n > 0); + + babl_polynomial_reset (poly, basis[0].scale); + + for (i = 0; i < basis_n; i++) + { + BablPolynomial temp; + + babl_polynomial_copy (&temp, &basis[i]); + babl_polynomial_scalar_mul (&temp, + babl_polynomial_gamma_inner_product (&temp, + gamma, + x0, x1)); + babl_polynomial_add (poly, &temp); + } +} + +static void +babl_polynomial_gram_schmidt (BablPolynomial *basis, + int basis_n, + double x0, + double x1) +{ + int i; + + for (i = 0; i < basis_n; i++) + { + if (i > 0) + { + BablPolynomial temp; + + babl_polynomial_project_copy (&temp, &basis[i], basis, i, x0, x1); + babl_polynomial_sub (&basis[i], &temp); + } + + babl_polynomial_normalize (&basis[i], x0, x1); + } +} + +static void +babl_polynomial_basis (BablPolynomial *basis, + int basis_n, + int scale) +{ + int i; + + for (i = 0; i < basis_n; i++) + { + babl_polynomial_set_degree (&basis[i], i, scale); + + basis[i].coeff[0] = 1.0; + memset (&basis[i].coeff[1], 0, i * sizeof (double)); + } +} + +static void +babl_polynomial_orthonormal_basis (BablPolynomial *basis, + int basis_n, + double x0, + double x1, + int scale) +{ + babl_polynomial_basis (basis, basis_n, scale); + babl_polynomial_gram_schmidt (basis, basis_n, x0, x1); +} + +void +babl_polynomial_approximate_gamma (BablPolynomial *poly, + double gamma, + double x0, + double x1, + int degree, + int scale) +{ + BablPolynomial *basis; + + babl_assert (poly != NULL); + babl_assert (gamma >= 0.0); + babl_assert (x0 < x1); + babl_assert (degree >= BABL_POLYNOMIAL_MIN_DEGREE && + degree <= BABL_POLYNOMIAL_MAX_DEGREE); + babl_assert (scale >= BABL_POLYNOMIAL_MIN_SCALE && + scale <= BABL_POLYNOMIAL_MAX_SCALE); + + basis = alloca ((degree + 1) * sizeof (BablPolynomial)); + + babl_polynomial_orthonormal_basis (basis, degree + 1, x0, x1, scale); + + babl_polynomial_gamma_project_copy (poly, gamma, basis, degree + 1, x0, x1); + babl_polynomial_shrink (poly); +} + +#endif diff --git a/babl/babl-polynomial.h b/babl/babl-polynomial.h new file mode 100644 index 0000000..d883542 --- /dev/null +++ b/babl/babl-polynomial.h @@ -0,0 +1,79 @@ +/* babl - dynamically extendable universal pixel conversion library. + * Copyright (C) 2017, Øyvind Kolås and others. + * + * babl-polynomial.h + * Copyright (C) 2017 Ell + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 3 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General + * Public License along with this library; if not, see + * . + */ + +#ifndef _BABL_POLYNOMIAL_H +#define _BABL_POLYNOMIAL_H + + +/* BablPolynomial is an opaque type representing a real polynomial of a real + * variable. In addition to a degree, polynomials have an associated *scale*, + * which divides the exponents of the polynomial's terms. For example, a + * polynomial of degree 3 and of scale 1 has the form + * `c0*x^0 + c1*x^1 + c2*x^2 + c3*x^3`, while a polynomial of degree 3 and of + * scale 2 has the form `c0*x^0 + c1*x^0.5 + c2*x^1 + c3*x^1.5`. + */ + + +#define BABL_POLYNOMIAL_MIN_DEGREE 0 +#define BABL_POLYNOMIAL_MAX_DEGREE 10 + +#define BABL_POLYNOMIAL_MIN_SCALE 1 +#define BABL_POLYNOMIAL_MAX_SCALE 2 + + +typedef struct BablPolynomial BablPolynomial; + +typedef double (* BablPolynomialEvalFunc) (const BablPolynomial *poly, + double x); + + +struct BablPolynomial +{ + BablPolynomialEvalFunc eval; + int degree; + int scale; + double coeff[BABL_POLYNOMIAL_MAX_DEGREE + 1]; +}; + + +/* evaluates `poly` at `x`. */ +static inline double +babl_polynomial_eval (const BablPolynomial *poly, + double x) +{ + return poly->eval (poly, x); +} + + +/* calculates the polynomial of maximal degree `degree` and of scale `scale`, + * that minimizes the total error w.r.t. the gamma function `gamma`, over the + * range `[x0, x1]`. + */ +void +babl_polynomial_approximate_gamma (BablPolynomial *poly, + double gamma, + double x0, + double x1, + int degree, + int scale); + + +#endif -- 2.30.2